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ABSTRACT 

This  paper  outlines  a  methodology  for  determining 
the  statistics  associated  with  the  time  evolution  of  a 
nonlinear  multi-body  dynamic  system  operated  under 
input  uncertainty.  The  focus  is  on  the  dynamics  of  ground 
vehicle  systems  in  environments  characterized  by  multiple 
sources  of  uncertainty:  road  topography,  friction 
coefficient  at  the  road/tire  interface  and  aerodynamic 
force  loading.  Drawing  on  parametric  maximum 
likelihood  estimation,  the  methodology  outlined  is  general 
and  can  be  applied  to  systematically  study  the  impact  of 
sources  of  uncertainty  characterized  herein  by  random 
processes.  The  proposed  framework  is  demonstrated 
through  a  study  that  characterizes  the  uncertainty  induced 
in  the  loading  of  the  lower  control  arm  of  an  SUV  type 
vehicle  by  uncertainty  associated  with  road  topography. 

1.  INTRODUCTION 

The  goal  of  this  work  is  to  establish  an  analytically 
sound  and  computationally  efficient  framework  for 
quantifying  uncertainty  in  the  dynamics  of  complex  multi¬ 
body  systems.  The  motivating  question  for  this  effort  is  as 
follows:  how  can  one  predict  an  average  response  and 
produce  a  confidence  interval  in  relation  to  the  time 
evolution  of  a  complex  multi-body  system  given  a  certain 
degree  of  input  uncertainty?  Herein,  of  interest  is 
answering  this  question  for  ground  vehicle  systems  whose 
dynamics  are  obtained  as  the  solution  of  a  set  of 
differential-algebraic  equations  (Hairer;  Wanner  1996). 
The  differential  equations  follow  from  Newton's  second 
law.  The  algebraic  equations  are  nonlinear  kinematic 
equations  that  constrain  the  evolution  of  the  bodies  that 
make  up  the  system  (Haug  1989). 

The  motivating  question  above  is  relevant  for  vehicle 
Condition-Based  Maintenance  (CBM)  where  the  goal  is  to 
predict  durability  and  fatigue  of  system  components.  The 
statistics  of  lower  control  arm  loading  in  a  High-Mobility 
Multi- Wheeled  Vehicle  (HMMWV)  obtained  through  a 
multi-body  dynamics  simulation  become  the  input  to  a 


durability  analysis  that  can  predict  in  a  stochastic 
framework  the  condition  of  the  part  and  recommend  or 
postpone  system  maintenance.  A  stochastic 
characterization  of  system  dynamics  is  also  of  interest  in 
understanding  the  limit  behavior  of  a  dynamic  system.  For 
instance,  providing  real-time  confidence  intervals  for 
certain  vehicle  maneuvers  are  useful  in  assessing  its 
control  when  operating  in  icy  road  conditions. 

Vehicle  dynamics  analysis  under  uncertain 
environment  conditions,  e.g.  road  profile  (elevation, 
roughness,  friction  coefficient)  and  aerodynamic  loading, 
requires  approaches  that  draw  on  random  functions.  The 
methodology  is  substantially  more  involved  than  required 
for  handling  uncertainty  that  enters  the  dynamic  response 
through  discrete  design  parameters  associated  with  the 
model.  For  instance,  uncertainty  in  suspension  spring 
stiffness  or  damping  rates  can  be  handled  through  random 
variables.  In  this  case,  methods  such  as  the  polynomial 
chaos  (PC),  see,  for  instance,  (Xiu;  Karniadakis  2002)  are 
suitable  provided  the  number  of  random  variables  is 
small.  This  approach  is  not  suitable  here  since  a 
discretization  of  the  road  leads  to  a  very  large  number  of 
random  variables  (the  road  attributes  at  each  road  grid 
point).  Moreover,  the  PC  methodology  requires  direct 
access  and  modification  of  the  computer  program  used  to 
run  the  deterministic  simulation  of  the  dynamic  system  to 
produce  first  and  second  order  moment  statistical 
information.  This  represents  a  serious  limitation  if  relying 
on  commercial  off-the-shelf  (COTS)  software,  which  is 
most  often  the  case  in  industry  when  running  complex 
high-fidelity  vehicle  dynamics  simulation. 

In  conjunction  with  Monte  Carlo  analysis,  the 
alternative  considered  herein  relies  on  random  functions 
to  capture  uncertainty  in  system  parameters  and/or  input. 
Limiting  the  discussion  to  three-dimensional  road  profiles, 
the  methodology  samples  a  posterior  distribution  that  is 
conditioned  on  available  road  profile  measurements.  Two 
paths  can  be  followed  to  implement  this  methodology. 
The  first  draws  on  a  parametric  representation  of  the 
uncertainty;  the  second  is  nonparametric  in  nature  and  as 
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such  is  more  expensive  to  implement.  It  can  rely  on 
smoothing  techniques  for  kernel  estimation  such  as 
Nadaray  a- Watson  estimation,  see,  for  instance, 
(Wasserman  2006),  or  draw  on  the  orthogonal 
discretization  of  the  spectral  representation  of  positive 
definite  functions  as  suggested,  for  instance,  by  Genton 
and  Gorsich  (2002).  The  parametric  approach  is  used  in 
this  paper  by  considering  Gaussian  Random  Functions  as 
priors  for  the  road  profiles.  Furthermore,  the  discussion 
will  be  limited  to  stationary  processes  although 
undergoing  research  is  also  investigating  the  nonstationary 
case. 

As  is  always  the  case,  the  use  of  a  parametric  model 
raises  two  legitimate  questions:  why  a  particular 
parametric  model,  and  why  is  it  fit  to  capture  the  statistics 
of  the  problem.  Gaussian  Random  Functions  (GRF)  are 
completely  defined  by  their  correlation  function,  also 
known  as  variogram  (Adler  1990;  Cramer;  Leadbetter 
1967).  Consequently,  scrutinizing  the  choice  of  a 
parametric  GRF  model  translates  into  scrutinizing  the 
choice  of  correlation  function.  There  are  several  families 
of  correlation  functions,  the  more  common  being 
exponential,  Matern,  linear,  spherical,  and  cubic  (see,  for 
instance  Santner  et  al.  (2003)).  In  this  context,  and  in 
order  to  demonstrate  the  proposed  framework  for 
uncertainty  quantification  in  multi-body  dynamics,  a 
representative  problem  will  be  investigated  in  conjunction 
with  the  selection  of  a  GRF-based  prior.  Specifically,  an 
analysis  will  be  carried  out  to  assess  the  sensitivity  of  the 
response  of  a  vehicle  to  uncertainty  in  system  input,  here  a 
road  profile.  Of  interest  is  the  load  history  for  the  lower- 
control  arm  of  an  HMMWV,  a  key  quantity  in  the  CBM 
of  the  vehicle.  The  parametric  priors  considered  are  (/)  a 
GRF  with  a  squared  exponential  correlation  function,  and 
(ii)  the  Ornstein-Uhlenbeck  process.  Pronounced 
sensitivity  of  the  statistics  of  the  loads  acting  on  the  lower 
control  arm  with  respect  to  the  choice  of  parametric  model 
would  suggest  that  serious  consideration  needs  to  be  given 
to  the  nonparametric  route,  where  the  empirical  step  of 
variogram  selection  is  avoided  at  the  price  of  a  more 
complex  method  and  increase  in  simulation  time. 

2.  UNCERTAINTY  HANDLING 
METHODOLOGY 

The  discussion  herein  concerns  handling  uncertainty 
in  spatial  data.  This  situation  commonly  arises  when 
limited  information  is  used  to  generate  road  profiles 
subsequently  used  in  the  dynamic  analysis  of  a  ground 
vehicle.  The  uncertainty  handling  in  aerodynamic  loads, 
which  can  be  addressed  similarly,  is  not  of  primary 
interest  in  this  study  and  will  be  omitted. 

The  uncertainty  quantification  framework  proposed  is 
described  in  Figure  1 .  The  assumption  is  that  learning  data 
is  available  as  the  result  of  field  measurements. 


Figure  1 :  Proposed  uncertainty  quantification  framework 

Referring  to  Figure  2,  the  measured  data  is  provided  on  a 
“coarse”  measurement  grid:  at  each  (xvx2)  location  an 
elevation  y  is  available.  From  an  uncertainty 
characterization  perspective,  the  dimension  of  this 
problem  is  d  =  2  .  For  dynamic  analysis,  road  information 
is  ideally  available  everywhere  on  the  road  as  a 
continuous  data.  As  this  is  not  possible,  data  is  provided 
on  a  fine  grid  (right  image  in  Figure  2).  If  working  with  a 
parametric  model,  a  correlation  function  is  selected  and  a 
learning  stage  follows.  Its  outcome,  a  set  of  hyper¬ 
parameters  associated  with  the  correlation  function,  is 
instrumental  in  generating  the  mean  and  covariance  matrix 
ready  to  be  used  to  generate  sample  road  surfaces  on  the 
user  specified  fine  grid. 


Figure  2.  Coarse  grid  for  learning,  and  fine  grid  employed 
in  sampling  for  Monte  Carlo  analysis.  Here  d -  2. 

Gaussian  Random  Functions  (GRF)  or  processes  lead 
to  a  very  versatile  approach  for  the  simulation  of  infinite 
dimensional  uncertainty.  By  definition,  a  spatially 
distributed  random  function  y(x) ,  xe  Rd ,  is  a  GRF  with 
mean  function  m(x;^)  and  correlation  function 
k(x,x;02)  for  any  set  of  space  points 
X  =  {xPx2,...,xM}, 
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y(X)=  ~W(m(X;01),K(X,X;02)).  (1) 

v  y  m  j 

Here  yj  =  y(x.) ,  me  Mm  ,  Kg  RMxM  ?  and  N( m,K)  is 

the  M-variate  normal  distribution  with  mean  m  and 
covariance  K  given  by 

m(xx,0x) 
m(x2,0x) 

m(xM,Ox) 

k(xx,x';02)  •••  k{xx,xN\02) 

k(x2,x'x;02)  •••  k(x2,xN\02) 

k(xM,x[\02)  •••  k(xM,xN\02) 

where  X'  =  {x',  x2 , . . . ,  }  .  The  hyper-parameters  0X  and 
02  associated  with  the  mean  and  covariance  functions  are 
obtained  from  a  data  set  y  (D)  at  nodes  D  =  { dx , . . . ,  dM  }  . 
The  posterior  distribution  of  the  variable  y (S)  at  node 
points  S  =  {sj , . . . ,  sN  }  ,  consistent  with  y  (D) ,  is 
N{ f*,K*)  (Rasmussen;  Williams  2006),  where 


K  (X,X';fl2)  = 


m(X;0,)  = 


Here  W  =  y  -m(0)  and  y  is  the  observed  data.  Note  that 
0  =  {0X,02}  >  and  ^e  dependence  on  the  hyper-parameters 

0  appears  by  means  of  the  coordinates  x  .  The  gradients 
of  the  likelihood  function  can  be  computed  analytically 
(Rasmussen;  Williams  2006): 

Y 

m(<9)  K(6>r1W 

y 

-  log  p(y  1 0)  =  1  tr  (K  ((9)-i  w(K  (6>)_1  W)r  -  K(6»)_1 ) 

dd2j  2  dd2j  J 

MATLAB’s  f solve  function,  which  implements  a  quasi- 

Newton  approach  for  nonlinear  equations,  was  used  to 

solve  the  first  order  optimality  conditions 

5logp(yl0)  _  ,  dlogp(y\0)  _  t  t 

- - - =  0  and  - - - =  0  and  determine 

do{j  de2j 

the  hyper-parameters  6X  and  02 .  The  entire  approach 
hinges  at  this  point  upon  the  selection  of  the  parametric 
mean  and  covariance  function.  It  is  common  to  select  a 
zero  mean  prior  m  =  0 ,  in  which  case  only  the  02j  hyper¬ 
parameters  associated  with  the  covariance  matrix  remain 
to  be  inferred  through  MLE. 

2.2  Covariance  Function  Selection 


f *  =  K(S,  D;  02  )K’*  (D,  D;  02 )  (y  (D)  -  m(£>;  d, ))  +m(,S;  dl ) 
K*  =  K(S,S\02)-K(S,D\62)K-\D,D-,02)K(D,S-,02) 

The  key  issues  in  sampling  from  this  posterior  are  a) 
how  to  obtain  the  hyper-parameters  from  data,  and  b)  how 
to  sample  from  A(f*,K*),  especially  in  the  case  where 
M  is  very  large.  The  classical  way  to  sample  relies  on  a 
Cholesky  factorization  of  K* ,  a  costly  order  0(M3) 
operation.  The  efficient  sampling  question  is  discussed  at 
length  in  (Schmitt  et  al.  2008a).  A  brief  description  of  the 
hyper-parameter  calculation  follows. 

2.1  Parameter  Estimation 


The  parametric  covariance  function  adopted 
determines  the  expression  of  the  matrix  K(0)  of  the 
previous  subsection,  and  it  requires  an  understanding  of 
the  underlying  statistics  associated  with  the  data.  In  what 
follows,  the  discussion  focuses  on  four  common  choices 
of  correlation  function:  squared  exponential  (SE), 
Ornstein-Uhlenbeck  (OU)  (Uhlenbeck;  Ornstein  1930), 
Matern  (Matern  1960),  and  neural  network  (NN)  (Neal 
1996). 

The  SE  correlation  function  assumes  the  form 
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2  /y 

(x2-x2) 

2/y\ 
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e2l 

^22 

y 

The  method  used  herein  for  the  estimation  of  the 
hyper-parameters  data  draws  on  maximum  likelihood 
estimation  (MLE)  (Rasmussen;  Williams  2006). 
Specifically,  it  selects  those  hyper-parameters  that 
maximize  the  log-likelihood  function  associated  with  the 
measured  data.  In  the  multivariate  Gaussian  with  mean 
m(<9)eRM  and  covariance  matrix  K(#)eRMxM  case, 
the  log-likelihood  function  assumes  the  form 

log  p( y  I  6)  =  WrK(6»)_1  W  -  hog  I  K(6»)  I  ~  log  2^ 


where  y  =  1 .  The  hyper-parameters  02X  and  022  are 
called  the  characteristic  lengths  associated  with  the 
stochastic  process  and  control  the  degree  of  spatial 
correlation.  Large  values  of  these  coefficients  lead  to  large 
correlation  lengths,  while  small  values  reduce  the  spatial 
correlation  leading  in  the  limit  to  white  noise,  that  is, 
completely  uncorrelated  data.  The  SE  is  the  only 
continuously  differentiable  member  of  the  family  of 
exponential  GRF.  As  such,  it  is  not  commonly  used  for 
capturing  road  profiles,  which  are  typically  not 
characterized  by  this  level  of  smoothness.  Rather,  Stein 
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(Stein  1999)  recommends  the  Matern  family  with  the 
correlation  function 


k(r\02) 


2l-v 

nv) 


s/2vA 

V 

Kv 

1  ) 

l  1  ) 

(5) 


where  Kv  is  the  modified  Bessel  function  and  v  and  / 
are  positive  hyper-parameters.  The  degree  of  smoothness 
of  the  ensuing  GRF  can  be  controlled  through  the 
parameter  v  :  the  GRF  is  p  -times  differentiable  if  and 
only  if  v  >  p  .  Note  that  selecting  y  =  2  in  Eq.  (4)  leads 
to  the  OU  random  process,  which  is  also  a  nonsmooth 
process  although  not  as  versatile  as  the  Matern  family. 

The  three  covariance  models  discussed  so  far:  SE, 
OU,  and  Matern  are  stationary.  Referring  to  Eq.  (1),  this 
means  that  for  any  set  of  points  X  =  {xpx2,...,xM}  , 

where  M  is  arbitrary,  and  for  any  vector  he  ,  y(X) 
and  y(X  +  h)  always  have  the  same  mean  and  covariance 
matrix.  In  particular,  when  M  - 1 ,  this  means  that  the 
GRF  should  have  the  same  mean  and  variance 
everywhere.  Clearly,  the  stationary  assumption  does  not 
hold  in  many  cases.  For  vehicle  simulation,  consider  the 
case  of  a  road  with  a  pothole  in  it,  which  cannot  be 
captured  by  stationary  processes.  A  versatile 
nonstationary  neural  network  covariance  function  has 
been  proposed  by  Neal  (1996): 

k(x,x';£)  =  —  sin 
n 

V  V  '  '  '  '  J 


2xr£x' 


J(  1  +  2xrZx)(l  +  2x/rZx/) 


(6) 


where  x  =  (l,  xl , . . . ,  xd  )T  is  an  augmented  input  vector;  the 

symmetric  positive  definite  matrix  L  contains  the 
parameters  associated  with  this  GRF  that  are  determined 
through  MLE.  Note  that  for  the  road  profile  problem 
d  =  2. 

Using  this  parameter  estimation  approach,  a  mean 
and  co-variance  function  for  Gaussian  processes  is 
determined.  This  is  then  used  to  generate  new  roads  which 
are  statistically  equivalent  to  the  road  used  in  the  learning 
process.  Figure  5.  shows  the  vehicle  model  in 
ADAMS/Car  on  two  such  roads. 


3.  NUMERICAL  EXPERIMENT:  PARAMETRIC 
MODEL  SENSITIVITY 

The  numerical  experiments  carried  out  illustrate  how 
the  proposed  uncertainty  quantification  framework  is  used 
to  predict  an  average  behavior  and  produce  a  confidence 
interval  in  relation  to  the  time  evolution  of  a  nonlinear 
multi-body  system.  A  high-fidelity  vehicle  model  is 
considered  and  its  time  evolution  is  marked  by  uncertainty 
stemming  from  measurements  of  the  road  profile.  This 
setup  was  chosen  due  to  its  relevance  in  CBM,  where  the 


interest  is  the  statistics  of  the  loads  acting  on  the  vehicle 
for  durability  analysis  purposes.  Note  that  a  similar 
analysis  is  carried  out  for  a  simplified  scenario  that  does 
not  involve  the  MLE  learning  stage  by  Schmitt  et  al. 
(2008b)  in  conjunction  with  quantifying  the  uncertainty  of 
vehicle  dynamics  when  running  on  icy  roads  with  a 
stochastic  distribution  of  the  tire/road  friction  coefficient. 

3.1  Vehicle  Model 

The  vehicle  of  interest  in  this  work  is  a  SUV-type 
vehicle  similar  to  the  Army’s  High  Mobility  Multi- 
Wheeled  Vehicle  (HMMWV).  A  high-fidelity  model  of 
the  vehicle  was  generated  in  ADAMS,  a  widely  used 
COTS  software  that  contains  a  template  library, 
ADAMS/Car,  dedicated  to  ground  vehicle  modeling  and 
simulation.  The  steering  system  is  of  rack-and-pinion 
type.  The  vehicle  is  equipped  with  an  Ackerman  type 
suspension  system.  The  front  and  rear  suspensions  have 
the  same  topology  but  different  link  lengths.  The  location 
of  the  suspension  subsystem  is  parameterized  with  respect 
to  the  chassis  of  the  vehicle  to  allow  for  an  easy  editing  of 
the  assembly  topology.  Although  not  a  topic  of  interest  in 
this  paper  that  deals  with  infinite  dimensional  stochastic 
processes,  this  parameterization  of  the  suspension  allows 
for  a  Monte-Carlo  based  approach  to  quantify  the 
uncertainty  in  vehicle  dynamics  produced  by  this  model 
subcomponent. 
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Figure  3.  Up:  Vehicle  model,  (no  chassis  geometry 
shown).  Down:  Schematic  of  how  different  subsystems  are 
assembled  to  make  a  full  vehicle  model. 
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Figure  3  shows  the  topology  of  a  vehicle  with  front 
and  rear  suspension,  wheels,  and  steering  subsystems. 
Different  vehicle  subsystems  are  individually  modeled  and 
are  integrated  together  to  form  a  full  vehicle  model.  The 
chassis  of  the  vehicle  is  modeled  as  a  single  component 
having  appropriate  mass-inertia  properties.  ADAMS/Flex 
is  typically  used  to  add  compliance  to  the  steering  and/or 
chassis  components  of  the  vehicle.  This  process  makes  use 
of  modal  neutral  format  (MNF)  files  created  using  a  third 
party  finite  element  package  such  as  ABAQUS. 

3.2  Tire  and  Road  Models 

Two  types  of  external  forces  act  under  normal 
conditions  on  the  vehicle  and  influence  its  dynamics: 
forces  at  the  tire-road  interface  and  aerodynamic  forces. 
Given  the  range  of  speeds  at  which  the  considered  ground 
vehicle  is  driven,  the  former  forces  are  prevalent  and  high- 
fidelity  simulation  depends  critically  on  their  accurate 
characterization.  In  this  context,  the  SUV  configuration  of 
interest  is  designed  to  drive  over  a  variety  of  terrains, 
from  flat  smooth  pavement  to  off-road  conditions. 

In  this  study,  the  tire  is  modeled  using  a  high-fidelity 
COTS  software  called  FTire  (Gipser  2005).  FTire  can  be 
used  for  vehicle  handling,  ride  comfort  and  durability 
studies  on  even  or  uneven  roadways  with  extremely  short 
obstacle  wavelengths.  It  is  a  physics-based,  highly 
nonlinear,  dynamic  tire  model,  which  is  fast  (typically 
only  5  to  20  times  slower  than  real-time)  and  numerically 
robust.  The  tire  belt  is  described  as  an  extensible  and 
flexible  ring  carrying  bending  loads,  elastically  founded 
on  the  rim  by  distributed,  partially  dynamic  stiffness  in  the 
radial,  tangential,  and  lateral  directions.  The  tire  model  is 
accurate  at  relatively  high  frequencies  (up  to  120  Hz)  both 
in  longitudinal  and  lateral  directions.  It  works  out  of,  and 
up  to,  a  complete  standstill  without  any  additional 
computing  effort  or  model  switching.  It  is  suitable  for 
demanding  applications  such  as  ABS  braking  on  uneven 
roadways. 

The  road  models  supported  should  accommodate  the 
fidelity  representation  level  required  by  the  tire  model 
considered.  The  road  modeling  environment  of  choice  is 
the  one  provided  by  ADAMS,  where  the  road  is  defined 
by  a  text  based  data  file  (rdf).  This  rdf  file  contains  the 
information  about  road  size,  type  (flat,  periodic  obstacles, 
stochastic  3D)  and  coefficients  of  friction  over  the  road 
surface.  Defining  a  flat  road  is  trivial  and  an  obstacle  in 
the  road  (curb,  roof-shaped)  can  be  defined  by  specifying 
the  size  and  shape  of  the  obstacle.  For  roads  with  varying 
elevations  in  both  lateral  and  longitudinal  directions,  a 
tessellated  road  definition  is  supported  in  ADAMS.  The 
tessellated  road  is  described  by  a  set  of  vertices/nodes, 
which  are  grouped  in  sets  of  three  to  create  a  triangulated 
mesh  that  describes  the  entire  road  surface  (see  Figure  4). 
A  coefficient  of  friction  can  be  specified  for  each  triangle. 


This  road  definition  works  well  with  both  synthesized  and 
measured  road  data  and  it  was  the  solution  embraced  for 
the  numerical  results  reported  in  this  work. 

The  disadvantage  of  this  road  definition  is  that  for 
each  time  step,  the  simulation  has  to  check  each  triangle 
for  contact  with  the  tire  patch  leading  to  a  major 
computational  bottleneck  (extremely  long  simulation 
times  for  large  road  profiles  with  high  resolution).  In 
order  to  simulate  large  road  profiles,  the  regular  grid  road 
(rgr)  file  format  can  be  used  in  conjunction  with  FTire.  A 
conversion  from  tessellated  to  the  new  rgr  format  leads  to 
smaller  file  sizes  and  significantly  reduces  CPU  time  per 
simulation  step  in  that  the  CPU  time  required  for  the 
tire/terrain  interaction  is  independent  of  the  dimensions 
(length/width)  of  the  road  profile.  This  approach  has  been 
used  in  the  past  to  analyze  ride  maneuvers  that  cover  mile 
long  distances. 


Figure  4:  Road  tessellation  through  a  triangulation 
applied  to  a  road  profile 

3.3  Numerical  Results,  Square  Exponential 

The  SUV  model  discussed  in  subsection  3.1  was 
equipped  with  a  set  of  four  tires  generated  in  the  FTire 
modeling  and  simulation  package  discussed  in  subsection 
3.2.  The  vehicle  model  was  exercised  through  a  straight- 
line  maneuver  over  a  road  profile  for  which  information  is 
available  on  a  grid  as  follows  (see  also  Figure  4):  in  the  x- 
direction,  information  is  provided  every  0.25  feet  in  180 
slices.  In  the  y-direction,  the  data  is  provided  at  a  distance 
of  four  feet  apart  in  three  slices.  The  length  of  the  course 
in  the  x-direction  was  approximately  45  feet.  The  width  of 
the  road  was  8  feet.  Although  not  reported  here, 
simulations  up  to  one  mile  long  have  been  run  using  this 
vehicle  configuration. 

The  stochastic  analysis  proceeded  according  to  the 
work  flow  in  Figure  1.  A  road  profile  was  provided  and 
considered  the  outcome  of  a  set  of  field  measurements  on 
a  180X3  grid  as  indicated  above.  MLE  was  carried  out; 
the  resulting  characteristic  lengths  were  Gx  =4.5355  and 
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Gy  =  0.8740  ,  which  are  identified  in  Eq.  (4)  with  02l  and 
022 ,  respectively.  A  set  of  200  road  profiles  were 
generated  by  sampling  of  the  posterior;  Figure  5  illustrates 
two  of  them.  The  road  profiles  were  relatively  smooth  in 
the  sense  that  there  was  no  road  geometric  feature  of 
length  comparable  to  the  length  of  the  tire/road  contact 
patch.  The  terrain  was  considered  rigid  and  with  a 
constant  friction  coefficient. 


Figure  5.  Vehicle  on  2  different  road  surfaces.  The  two 
roads  are  generated  using  the  same  training  data. 
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Figure  6.  Vehicle  response  obtained  using  ADAMS.  A 
subset  of  10  out  of  the  200  simulations  used  to  generate 
the  statistics  of  the  normal  force  response  is  displayed. 
The  plot  shows  the  reaction  force  in  a  suspension  joint. 

A  batch  of  200  ADAMS  simulations  were 

subsequently  carried  out  to  determine  the  statistics  of  the 
vehicle  response.  Ten  responses  are  illustrated  in  Figure  6, 
which  reports  the  reaction  force  in  a  suspension  joint.  Of 
interest  here  is  the  loading  in  the  lower  control  arm  (FCA) 
of  the  vehicle  suspension,  which  is  measured  by  the  loads 
experienced  by  the  suspension  bushings  connecting  the 
FCA  and  the  chassis.  There  are  two  such  bushings,  and 
there  are  two  more  connecting  the  upper  control  arm 
(UCA)  to  the  chassis  for  a  total  of  16  bushing  elements 
(eight  UCA  and  eight  FCA).  Average  behavior  and  a  95% 


confidence  interval  are  provided  for  the  FCA  bushing 
load  in  Figure  7.  Note  that  all  results  reported  in  the  plots 
are  in  SI  units. 


Figure  7.  Statistics  of  vertical  load,  bushing  attached  to 
FCA. 


Figure  8.  Statistics  of  vertical  load,  bushing  attached  to 
FCA.  Detailed  view  focused  on  the  last  part  of  run. 


Figure  9.  Statistics  of  vertical  acceleration  as  measured  at 
the  CG  position  of  the  chassis:  mean  and  95%  confidence 
interval. 
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3.4  Numerical  Results,  Ornstein-Uhlenbeck 

The  results  reported  in  this  subsection  are  obtained  by 
setting  y  =  2  in  Eq.  (4).  This  change  leads  to  Ornstein- 
Uhlenbeck  GRF  posteriors  that  lack  differentiability.  In 
fact,  the  SE  ( y  =  1 )  is  the  only  exponential  GRF  that  is 
continuously  differentiable.  Figure  10  shows  the  load 
history  for  the  same  bushing  element  that  was  considered 
for  the  results  in  Figure  7.  Figure  1 1  is  a  zoom-in  to  better 
gauge  the  95%  confidence  interval  for  the  FCA  load. 
Finally,  the  vertical  acceleration  associated  with  the  OU 
process  is  reported  in  Figure  12.  Note  that  for  OU  the  two 
characteristic  lengths  obtained  at  the  end  of  the  MFE 
stage  are  0x  =  1.0849e  +  003  and  6y  =0.1830  . 


Figure  10.  Vertical  load  statistics,  bushing  attached  to 
FCA. 


Figure  11.  Statistics  of  vertical  load,  bushing  attached  to 
FCA.  Detailed  view  focused  on  the  last  part  of  run. 

3.5  Discussion  of  Numerical  Results 

The  results  reported  in  Figure  8  and  Figure  1 1  suggest 
that  the  experimental  data  is  gathered  on  a  dense  enough 
grid.  Specifically,  there  is  relatively  small  variance  in  the 
response  of  the  vehicle,  which  is  a  very  desirable  response 
characteristic.  This  can  be  also  seen  in  Figure  8,  which 
illustrates  the  statistics  associated  with  the  last  part  of  the 
simulation.  This  is  an  indication  that  the  grid  used  is 
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Figure  12.  Statistics  of  vertical  acceleration,  as  measured 
at  the  CG  position  of  the  chassis:  mean  and  95% 
confidence  interval. 

sufficiently  fine  to  limit  the  amount  of  uncertainty 
stemming  from  the  uncertainty  in  road  profile.  In  other 
words,  the  measured  road  data  is  provided  at  a  level  of 
granularity  sufficient  to  pinpoint  with  good  precision  the 
load  history  for  the  force  acting  at  a  hot  point  of  the  FCA. 

The  results  reported  in  Figure  9  indicate  the  statistics 
associated  with  the  vertical  acceleration  of  the  center  of 
gravity  (CG)  of  the  chassis.  In  this  context,  information 
regarding  the  vertical  acceleration  and  the  jerk  measured 
at  a  location  where  the  vehicle  driver  is  positioned  is 
valuable  as  it  is  used  to  gauge  ride  comfort  and  the 
potential  for  vibration  induced  fatigue  during  long  term 
exposure  of  a  vehicle  driver. 

Finally,  the  comparison  of  results  reported  in  Figure  7 
and  Figure  10,  or  Figure  8  and  Figure  11,  or  Figure  9  and 
Figure  12,  demonstrates  the  qualitative  difference  between 
the  SE  and  OU  exponential  GRFs.  The  OU  family  leads  to 
processes  that  are  nonsmooth,  while  the  SE  family  leads 
to  road  profiles  that  are  continuously  differentiable.  This 
is  eventually  reflected  in  the  smoothness  of  the  output:  the 
SE  response  is  smoother  compared  to  the  OU  outcome. 
Nonetheless,  there  is  good  agreement  between  the  results 
obtained  with  the  OU  and  SE,  and  if  this  roughness 
associated  with  OU  leads  to  any  significant  change  in 
CBM  decisions  remains  to  be  investigated.  This  task  falls 
outside  the  scope  of  this  study. 

4.  CONCLUSIONS  AND  FUTURE  WORK 

This  paper  outlines  a  methodology  for  determining 
the  statistics  associated  with  the  time  evolution  of  a 
nonlinear  multi-body  dynamic  system  operated  under 
input  uncertainty.  The  focus  is  on  the  dynamics  of  ground 
vehicle  systems  in  environments  characterized  by  multiple 
sources  of  uncertainty:  road  topography,  friction 
coefficient  at  the  road/tire  interface  and  aerodynamic 
force  loading.  The  methodology  outlined  is  general  and 
can  be  applied  to  systematically  study  the  impact  of 
sources  of  uncertainty  that  were  characterized  herein  by 
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random  processes.  The  scope  of  the  discussion  is  limited 
to  the  case  of  unknown  road  profiles  at  the  wheel/road 
interface  of  an  SUV.  The  same  approach  can  be  used  in 
conjunction  with  vehicle  design  parameters  by  substituting 
the  GRF  machinery  with  that  of  sampling  from  a  random 
variable  distribution.  The  latter  scenario  is  simpler  and  not 
discussed  here.  It  would  closely  follow  the  methodology 
outlined  in  Figure  1  in  the  sense  that  one  could  deal  with 
parametric  or  nonparametric  models,  and  then  invoke 
MLE  to  obtain  the  hyper-parameters  associated  with  the 
posterior  distribution.  The  latter  is  subsequently  sampled 
in  a  Monte-Carlo  analysis  to  gauge  the  impact  of 
uncertainty  on  the  dynamics  of  the  model. 

The  results  reported  herein  suggest  that  the  choice  of 
correlation  function  is  important.  Further  work  is  needed 
to  better  understand  the  sensitivity  of  the  system  response 
with  respect  to  the  correlation  function.  In  this  context, 
two  directions  of  future  work  could  prove  insightful.  First, 
it  would  be  useful  to  bring  into  the  picture  other  GRF 
correlation  functions,  both  stationary  and  non-stationary, 
to  understand  the  importance  of  the  parametric  model 
choice.  This  direction  is  currently  under  investigation 
(Datar  2008).  Second,  it  would  be  useful  to  lift  the 
requirement  that  the  stochastic  approach  be  molded  upon 
the  GRF  idea.  Nonparametric  models  should  be 
considered  and  the  improved  flexibility  should  be 
weighted  against  their  more  involved  analytical 
derivation,  method  implementation,  and  computational 
burdens. 

Even  before  considering  nonparametric  techniques  the 
computational  load  associated  with  the  methodology 
proposed  can  be  significant.  This  is  particularly  the  case 
when  the  amount  of  learning  data  is  vast  and/or  the 
posterior  is  sampled  on  a  very  fine  grid.  This  issue  can  be 
addressed  by  either  allocating  more  computational  power 
to  solve  the  problem  or  by  considering  new  sampling 
techniques.  Recent  inroads  into  better  sampling 
techniques  have  been  recently  reported  by  Schmitt, 
Anitescu  et  al.  (2008a),  where  a  substantial  reduction  in 
sampling  effort  is  demonstrated  by  the  use  of  two  new 
approaches.  The  first  makes  a  periodicity  assumption  that 
enables  a  Fast  Fourier  Transform  (FFT)  technique  to  be 
used  for  posterior  sampling,  while  the  second  uses 
compact  kernel  covariance  functions  to  sample  the 
posterior  on  small  sliding  windows  that  are  continuously 
moved  to  follow  the  vehicle  in  its  motion  over  the  region 
of  interest. 
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